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Abstract 



u 

, Two approximations, derived from continuous expansions of Riemann-Liouville frac- 

tional derivatives into series involving integer order derivatives, are studied. Using those 
series, one can formally transform any problem that contains fractional derivatives into 
a classical problem in which only derivatives of integer order are present. Correspond- 
ing approximations provide useful numerical tools to compute fractional derivatives of 
functions. Application of such approximations to fractional differential equations and 
fractional problems of the calculus of variations are discussed. Illustrative examples show 
the advantages and disadvantages of each approximation. 
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1 Introduction 

Fractional calculus is the study of integrals and derivatives of arbitrary real or complex order. 
^ " Although the origin of fractional calculus goes back to the end of the seventeenth century, 

the main contributions have been made during the last few decades [311 132j . Namely it has 
proven to be a useful tool when applied to engineering and optimal control problems (see, e.g., 
|12 1 119 1 129])- There are several different definitions of fractional derivatives in the literature, 
such as Griinwald-Letnikov, Caputo, etc. Here we consider RiemanmLiouville fractional 
derivatives. 

Definition 1.1 (cf. |16|). Letx(-) be an absolutely continuous function in [a, b] andO < a < 1. 
Then, 



the left Riemann-Liouville fractional derivative of order a, a Df , is given by 

aD?x{t)= T(^)jtiy~ T) ~ ax{T)dT ' teM]; (i) 



'This is a preprint of a paper whose final and definite form will be published in: Asian Journal of Control. 
Submitted 13-Oct-2011; revised ll-Apr-2012; accepted 10-Aug-2012. 
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• the right Riemann-Liouville fractional derivative of order a, tD^, is given by 



,D?x(f) 



Due to the growing number of applications of fractional calculus in science and engineering 
(see, e.g., [SJ [33]), numerical methods are being developed to provide tools for solving 
such problems. Using the Griinwald-Letnikov approach, it is convenient to approximate the 
fractional differentiation operator, D a , by generalized finite differences. In [23] some problems 
have been solved by this approximation. In [10] a predictor-corrector method is presented 
that converts an initial value problem into an equivalent Volterra integral equation, while 
|17] shows the use of numerical methods to solve such integral equations. A good survey on 
numerical methods for fractional differential equations can be found in [13j . 

A new numerical scheme to solve fractional differential equations has been recently in- 
troduced in [7] and [15], making an adaptation to cover fractional optimal control problems. 
The scheme is based on an expansion formula for the Riemann-Liouville fractional derivative. 
Here we introduce a generalized version of that expansion and, together with a different ex- 
pansion formula that has been used to approximate the fractional Euler-Lagrange equation 
in [5], we perform an investigation of the advantages and disadvantages of approximating 
fractional derivatives by these expansions. The approximations transform fractional deriva- 
tives into finite sums containing only derivatives of integer order. We show the efficiency 
of such approximations to evaluate fractional derivatives of a given function in closed form. 
Moreover, we discuss the possibility of evaluating fractional derivatives of discrete tabular 
data. The application to fractional differential equations and the calculus of variations is 
also developed through some concrete examples. In each case we try to analyze problems 
for which the analytic solution is available. This approach gives us the ability of measuring 
the accuracy of each method. To this end, we need to measure how close we get to exact 
solutions. We use the 2-norm and define the error function E[x(-),x(-)] by 



where x(-) is defined on [a, b}. The results of the paper give interesting numerical procedures 
when applied to fractional problems of the calculus of variations. 

2 Expansion formulas to approximate fractional derivatives 

In this section two approximations for the left Riemann-Liouville derivative are presented. 
Both approximate the fractional derivatives by finite sums including only derivatives of integer 
order and are based on continuous expansions for the left Riemann-Liouville derivative. 

2.1 Approximation by a sum of integer order derivatives 

The right-hand side of (pQ) is expandable in a power series involving integer order derivatives 
[28]. Let (c,d), — oo < c < d < +oo, be an open interval in 1R, and [a, b] C (c, d) be such 
that for each t € [a, b] the closed ball Bb_ a (t), with center at t and radius b — o, lies in (c, d). 
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For any real analytic function x(-) in (c, d) we can give the following expansion formula: 

D"x(t) ~ V ("} {t ~ a)n ' a X M(t) where - - a) 

a D t x{t) - 2_, ^ r(n + 1 _ a) x it), where - ^ _ ^ + J} . (2) 

The condition Bb_ a (t) C (c, d) comes from the Taylor expansion of x(t — r) at t, for r G (a, t) 
and t G (a, 6). The proof of this statement that can be found in [28] uses a similar expansion 
for fractional integrals. Here we outline a direct proof due to our requirements in Section T2.31 
Since x(t) is analytic, it can be expanded as a convergent power series, i.e., 

' re! 

n=0 

and then by ([TJ 

^ {t)= w^nL ( (t - T)_<, S^"^ (f " T) T n t6<<i,,,) - (3) 

Termwise integration, followed by differentiation and simplification, leads to 

» D «*> = - + f(Tb) | ( (JX-D! + ^ xt " Ht)(t - 

and finally to expansion formula ([5]). From the computational point of view, one can take 
only a finite number of terms in ([2]) and use the approximation 

N / \ 1 

a A°x(t)^^C(re,a)(t-ar- a x(")(t), where C(n, a) = (" — — -. (4) 

' \re/l(re + l — a) 

n=o \ / v 1 / 

2.2 Approximation using moments of a function 

The following lemma gives the departure point to another expansion. For a proof see [9]. 

Lemma 2.1 (Lemma 2.12 of [9j). Let x(-) G AC [a, 6] and < a < 1. Then the left Riemann- 
Liouville fractional derivative a Dfx{-) exists almost everywhere in [a,b]. Moreover, a Dfx{-) G 
L p [a, b] for 1 < p < i and 



a D t a x(t) 



T(l-a) 



x(a) 



a ) J a 



te(a,b). (5) 



Let Vp(x(-)), p G N, denote the (p — 2)th moment of a function x(-) G ^4C 2 [a, 6] (cf. [7]): 
: = = (1 -p) [ (r - a) p - 2 x(r)dr, p G N, f > a. (6) 



Following [7J, it is easy to show that, by successive integrating by parts, (J5J> is reduced to 

^™ / \ x(a) , . „ i(a) . , r „ (t — a) 1_Q /"* / r — a\ 1_a . . , 

oA^W = rn J * " « + ^ i - a 1_Q + Vro S / 1 - 7 x(r)dT. 

r(l-a) r(2-a) r(2-a) J a \ t-a/ 

(7) 
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Using the binomial theorem we conclude that 
a D?x{t) 



x(a) , . „ It — a) 1 a . . , 
(t - a) a + E,„ ; x(a) 



Til -a) 



r(2 - a) 



(/-a) 1 "" [* (j^T(p-l + a) (t-o> p 



r(2-a) y a ^ r(a-i)p! 

Further integration by parts and simplification in (JSj) gives 



x(r)dT, t> a. (8) 



a £> t a x(t) = A(a)(t - a)" Q x(t) + B(a)(i - a)^ a i:(t) - ^ C(a,p)(t - af-^Vpit) , (9) 

p=2 



where V p (t) is defined by ([6]) and 

A(a) = 

= 

C(a,p) = 



1 



r(l-a) 

1 

r(2 - a) 



p=2 



r(p - 1 + a) 
r(a)(p-l)! 

r(p - 1 + a) 



1 



^ r(a-l)p! 
r(p - 1 + a) 



r(2-a)r(a-l) (p-1)! 



(10) 



The moments Vp(i), p = 2,3, . . ., are regarded as the solutions to the following system of 
differential equations: 

' v p (t) = (i- P )(t- a y~ 2 x(t) 

V p (a) = 0, p = 2,3,... 

For numerical purposes, only a finite number of terms in the series ([9|) are used. We approx- 
imate the fractional derivative as 

N 

a Dfx(t) ~A{a,N)(t-a)- a x(t) + B{a,N)(t-a) 1 - a x(t)-J2C(a,p)(t-a) 1 - p - a V p {t), (11) 

p=2 



where A(a,N) and B(a,N) are given by 
A(a,N) ' 



B(a,N) 



r(i-a) 
1 

r(2 - a) 



AT 

p=2 
N 



r(p - 1 + a) 
r(a)(p-l)! 

r(p - 1 + a) 



\ - i \p - i -r u 
^ r(a-l)p! 



P =i 



(12) 



(13) 



Remark 2.2. Our approximation (|lip is different from the one presented in J?]/: since the 
infinite series Y^Li T r(a-i)p\ tends to —1, B(a) = and thus 



N 



Q D?x(t) ~ ^l(a,iV)t- Q x(t) - J^C(a,p)* 1 - 1, - a ^(*). 

p=2 



(14) 
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However, regarding the fact that we use a finite sum, in practice one has 

N 



=1 



r(p - 1 + a) 
r(a - l)p! 



^0. 



Therefore, and similarly to \27j , we keep here the approximation in the form (|lip . TTie 
vaZue of B(a,N) for some values of N and for different choices of a is given in Tabled It 



N 


4 


7 


15 


30 


70 


120 


170 


5(0.1, N) 


0.0310 


0.0188 


0.0095 


0.0051 


0.0024 


0.0015 


0.0011 


5(0.3, N) 


0.1357 


0.0928 


0.0549 


0.0339 


0.0188 


0.0129 


0.0101 


5(0.5, N) 


0.3085 


0.2364 


0.1630 


0.1157 


0.0760 


0.0581 


0.0488 


5(0.7, iV) 


0.5519 


0.4717 


0.3783 


0.3083 


0.2396 


0.2040 


0.1838 


5(0.9, N) 


0.8470 


0.8046 


0.7481 


0.6990 


0.6428 


0.6092 


0.5884 


5(0.99, iV) 


0.9849 


0.9799 


0.9728 


0.9662 


0.9582 


0.9531 


0.9498 



Table 1: B(a,N) for different values of a and N. 

shows that even for a large N, when a tends to one, B(a,N) cannot be ignored. In Figured 
we plot B(a,N) as a function of N for different values of a. In Section [3] we compare both 



a=0.3 

0*0.5 
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ct=0.95 

ci.0.99 




100 120 140 160 180 



Figure 1: B{a,N) for different values of a and N. 

approximations with some examples. We also refer to f()\, [26^ . where such type of expansion 
formulas are studied. 

A similar argument gives the expansion formula for t-D" , the right Riemann-Liouville 
fractional derivative. We propose the following approximation: 

N 

t D%x(t) ~ A(a, N)(b - t)~ a x(t) - B(a, N)(b - tf- a x(t) - ^ C(a,p)(b - i) 1 " p ~ a W r p(<) , 

p=2 

where W p (t) = (1 - p) f^b - T) p ~ 2 x(r)dT. Here A(a,N) and B(a,N) are the same as CG]) 
and (fl~3|h respectively. 

Formula ([9]) consists of two parts: an infinite series and two terms including the first 
derivative and the function itself. It can be generalized to contain derivatives of higher-order. 
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Theorem 2.3. Fix n G N and let x(-) G C n [a,b]. Then, 



n-l 



a D?x(t) = — At - a)- a x{t) + V A(a,i)(t - a)^^«(t) 

r(i-a) ^ 



p=n 



-T(p — n + 1 + a) 



r(-a)r(l + a)(p-n + l)! 



(t - a)- Q x(t) + B(a,p)(t - a) n - l - p ~ a V p (t) 



, (15) 



where 



A(a, i) 
B(a,p) 



1 



r(* + 1 - a 



i+ E 



p=n— % 

F(p — n + 1 + a) 
r(-a)r(l + a)(p-ra + l)!' 



r(p — n + 1 + a) 
r(a — i)(p — n + i + 1)! 



, i = l,...,n-l, 



y p (t) = (p-n + 1) / (r - a) p - n x(r)dr. 



Proof. Successive integrating by parts in ([5]) gives 
aJ Dfx(t) 



" (a) (t - a)- + ; (a) ,(t - a) 1 - + • • • + ^""^^(t - a)- 1 - 



r(i - a 
1 



T(2 - a) 
(t-Tr^Wfr)*. 



r(n — a) 



r(i - a) J a 

Using the binomial theorem, we expand the integral term as 
r-t 



(t-r) 



n— I— a (n 



K ' ^ ' j^T{l-n + a)p\{t-a)r> J a V ; y ' 



Splitting the sum into p = and p = 1 . . . oo, and integrating by parts the last integral, we 
get 



™ / \ {t-a)~ a . . (t-a) n - 2 - a 

j£)f x(t) = ^ ; ; x(a) + • • • + v 



r(l-a) 

(t - a)"- 1 -" 



r(n — 1 — a) 



x^ n - 2 \a) 



T(n — a) 
T(n -l-a)2-t 



(n-2) (t) 



p=i x y 

r(p — n + 1 + a) 
, r(-n + 2 + a)(p-l)!(t-a)P J a 



{t -aY- l x^ n - l \T)dT. 



The rest of the proof follows a similar routine, i.e., by splitting the sum into two parts, the 
first term and the rest, and integrating by parts the last integral until x(-) appears in the 
integrand. □ 

Remark 2.4. The series that appear in A(a, i) are convergent for all i G {1, . . . , n — 1}. Fix 
an i and observe that 



E 



r(p - n + 1 + a) 



E 



T(p + a-i) 



— ' T(a — i)(p — n + i + 1)\ ^— ' r(a — 

n— r p=l 



L F (a -i, 1) - 1. 
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Since i > a, iFq(oc — i, 1) converges by Theorem 2.1.1 of [2]. In "practice we only use finite 
sums and for A(a, i) we can easily compute the truncation error. Although this is a partial 
error, it gives a good intuition of why this approximation works well. Using the fact that 
l-Fo(M) = i/a < (cf. Eq. (2.1.6) in we have 

1 T(p — n + 1 + a) 

r(i + 1 - a) p= ^ +1 T(a-i)(p-n + i + 1)1 

1 i z? i ■ -w v^ +1 r(p + a - in 
= r (j + i-a) E r (a - W j 

~ r(i + l-a) ^ r(a-i)p! 

In Ta&Ze we give some values for this error, with a = 0.5 and different values for i and 
N -n. 



i 





5 


10 


15 


20 


1 


-0.4231 


-0.2364 


-0.1819 


-0.1533 


-0.1350 


2 


0.04702 


0.009849 


0.004663 


0.002838 


0.001956 


3 


-0.007052 


-0.0006566 


-0.0001999 


-0.00008963 


-0.00004890 


4 


0.001007 


0.00004690 


0.000009517 


0.000003201 


0.000001397 



Table 2: Truncation errors of A(a, i, N) for a = 0.5. 



Remark 2.5. Using Euler's reflection formula, one can define B(a,p) of Theorem \2.3\ as 

-sin{-Ka)T{p — n + 1 + a) 



B(a,p) 



ir(p — n + 1) 



For numerical purposes, only finite sums are taken to approximate fractional derivatives. 
Therefore, for a fixed n £ N and N > n, one has 



n-l 



/V 



a D?x{t) Ri^2A(a,i,N)(t - af^x®^) + ^B{a,p){t - a) n - l - p - a V p (t), (16) 



i=0 



p=n 



where 



A(a,i,N) 
B(a,p) 



r(i + 1 - a) 



N 



i+E 



p=2 



T(p — n + 1 + a) 
T(a — i)(p — n + i + 1)! 



i = 0, . . . , n — 1, 



r(p — n + 1 + a) 



r(-a)r(l + a)(p-n + l) 



y p (t) = (p-n + 1) / (r - a) p - n x(r)dr. 



Similarly, we can deduce an expansion formula for the right fractional derivative. 
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Theorem 2.6. Fix n G N and x(-) G C n [a, 6]. T/ien, 



1 n— 1 

D?x{t) = T(l-a) {b ~ tr ° x(t) + 5 " *) 



*-<*_(*) 



E 



p—n 



i=l 

—T{p — n + 1 + a) 
T(-a)r(l + a)(p-n + 



— (6 - + B(a,p)(6 - tr^-PWpit) 



where 

A(a, i) 

B(a,p) 



-If 



T(i + l-a) 



» + E 



p=n— i 

(-l) n T(p-n + l + a) 
r(-a)r(l + a)(p-n + l)!' 

(p-n + 1) / (6 - T) p - n x(T)dT. 



r(p — n + 1 + a) 
r(— z + a)(p — n + 1 + i) 



, i = 1,. . . ,n- 1, 



W p (t) 

Proof. Analogous to the proof of Theorem 
2.3 Error estimation 



□ 



This section is devoted to the study of the error caused by choosing a finite number of terms 
in the expansions. For the expansion (0), we separate the error term in © and rewrite it as 

i d r>f {t _ Tra £(-irf(t) {t _ T)n ) iT 



aD?x(t) 



T(l — a) dt J a 



n=0 



J7! 



r( 



I^yl jf ((*-)- E 

' Ja \ n=N+l 



<-i)-^w (t _ r) .U (17) 



n! 



The first term in (|17p gives @ directly and the second term is the error caused by truncation. 
The next step is to give a local upper bound for this error, Et r (t). The series 

-l) n x^(t) 



E 

n=N+l 



-(t - r) n , r G (a,t), t G (a,b), 



is the remainder of the Taylor expansion of x(r) and thus bounded by 



(N+l)\ K l T > 



111 



which M = max 

T€[a,t] 



E tr (t) < 



M 



(r) 



. Then, 
d rt 



T{\- a){N + l)\dt 



(t-r) 



N+l-a 



dr 



M 



T(l-a)(N + l) 



-(t-a) 



N+l-a 



For approximation (jllh . we observe that the integrand in ((7|) can be expanded, by the 
binomial theorem, as 



1 



t — a 
t — a 



l-a 



r(p — 1 + a) ( t — a 



Ei^-it« 
„ T(a - Tjp\ \t-a 



p=0 

N 



^ Tjp-l + a) f r-a \ p 

g r( a-i)p! [t^) +Rn ^ 



(18) 



8 



where 



B M r(p-l+g) (r-ay 



p=N+l 



Substituting (fTKj) into (|7|), we get 

x(a) 



D?x(t) 



r(i-a 

+ r(2 - a) J„ 
x(a) 



( t _ Q )-a + ( f _ a )l-a 



N 

£ 

v p=0 



r(2 - a) 

r(p — 1 + a) ( t — a 



r(a — \ i — a 



+ i?Ar(r) i(r)dr 



r(i-a) 

(t-a) 1 - 
+ r(2 - a) /, 

(t-a) 1 " 



(t - a)~ a + * (o) v (t - a) 1 '" 



N 

£ 

v p=0 



r(2 - a) 

r(p — 1 + a) ( t — a 



r(a — \ t — a 



x(r)dr 



Rn{t)x(t)(It. 



T(2-a) 

At this point, we apply the techniques of [7] to the first three terms with finite sums. Then, 
we receive (HID with an extra term of truncation error: 



E tr (t) 



(t-a) 



l-a ft 



T(2-a) J a 
Since < < 1 for r £ [a, t], one has 

r(p - 1 + a 



Rn(t)x(t)cIt. 



\R n (t) 



s £ 

p=7V+l 

coo (1— o) 2 +l— Q 



< 



=7V F 



r(a - l)j>! 

dp 



£ 

p=7V+l 
p (l-a) 2 + l-a 



1 — a 
p 



s £ 

p=AT + l 



3 (l-a) 2 +l-a 



,2-Q 



(l-a)iV 



l-a ' 



Finally, assuming L2 = max |x(t)|, we conclude that 

r£[a,t] 



\E tr (t)\ <L 2 - 



,(l-a) 2 +l-a 



-(t-a) 



2-a 



'T(2-a){l-a)N 1 - 

In the general case, the error is given by the following result. 

Theorem 2.7. If we approximate the left Riemann-Liouville fractional derivative by the finite 
sum (|16p . then the error Et r (-) is bounded by 

e (n—l—a) 2 +n—l—a 

\ E tr(t)\ < L n TT7Z — - - a)»-« (19) 



where 



T(n - a)(n - 1 - a)N 
L, 



max 

T£[a,t] 



X {n) (T, 



From (|19p we see that if the test function grows very fast or the point t is far from a, 
then the value of N should also increase in order to have a good approximation. Clearly, if 
we increase the value of n, then we need also to increase the value of ./V to control the error. 
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3 Numerical evaluation of fractional derivatives 



In [23] a numerical method to evaluate fractional derivatives is given based on the Griinwald- 
Letnikov definition of fractional derivatives. It uses the fact that for a large class of functions, 
the Riemann-Liouville and the Griinwald-Letnikov definitions are equivalent. We claim that 
the approximations discussed so far provide a good tool to compute numerically the frac- 
tional derivatives of given functions. For functions whose higher-order derivatives are easily 
available, we can freely choose between approximations or (jlip . But in the case that dif- 
ficulties arise in computing higher-order derivatives, we choose the approximation (jlip that 
needs only the values of the first derivative and function itself. Even if the first derivative is 
not easily computable, we can use the approximation given by (|14p with large values for N 
and a not so close to one. As an example, we compute oDfx(t), with a = |, for x(t) = t A 
and x(t) = e 2t . The exact formulas of the derivatives are derived from 

oA°- 5 (n = r(n r +l + -o.5) tW " a5 and oA °' 5(eA4) = *-°- 5 ^-°.5(At), 

where E a a is the two parameter Mittag-Leffler function [24J. Figure [2] shows the results using 
approximation @. As we can see, the third approximations are reasonably accurate for both 
cases. Indeed, for x(t) = i 4 , the approximation with N = 4 coincides with the exact solution 
because the derivatives of order five and more vanish. The same computations are carried out 




(a) oA°- B (* 4 ) (b) oA°' 6 (e 2t ) 

Figure 2: Analytic (solid line) versus numerical approximation (|4"|) . 

using approximation (jlip . In this case, given a function x(-), we can compute V p by definition 
or integrate the system (jlOh analytically or by any numerical integrator. As it is clear from 
Figure [3l one can get better results by using larger values of N. Comparing Figures [2] and 
[3l we find out that the approximation @ shows a faster convergence. Observe that both 
functions are analytic and it is easy to compute higher-order derivatives. The approximation 
(jU) fails for non-analytic functions as stated in [7]. 

Remark 3.1. A closer look to Q and (jlip reveals that in both cases the approximations are 
not computable at a and b for the left and right fractional derivatives, respectively. At these 
points we assume that it is possible to extend them continuously to the closed interval [a, b] . 
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(a) oA°' 5 (t 4 ) (b) D°\e 2t ) 



Figure 3: Analytic (solid line) versus numerical approximation (lllj) . 

In what follows, we show that by omitting the first derivative from the expansion, as done 
in [7], one may loose a considerable accuracy in computation. Once again, we compute the 
fractional derivatives of x(t) = t 4 and x(t) = e 2 *, but this time we use the approximation given 
by (I14p . Figure H] summarizes the results. Our expansion gives a more realistic approximation 
using quite small N, 3 in this case. To show how the appearance of higher-order derivatives 
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Figure 4: Comparison of approximation (jllh proposed here and approximation (114D of [7] 

in generalization <\15h gives better results, we evaluate fractional derivatives of x(t) = i 4 and 
x(t) = e 2t for different values of n. We consider n = 1, 2, 3, N = 6 for x(t) = t 4 (Figure 5(a) ) 
and = 4 for x(t) = e 2t (Figure [5(b)] ) . 

3.1 Fractional derivatives of tabular data 

In many applications (see Section [5]) , the function itself is not accessible in a closed form, but 
as a tabular data for discrete values of the independent variable. Thus, we cannot use the 
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Figure 5: Analytic (solid line) versus numerical approximation (I15j) . 

definition to compute the fractional derivative directly. Our approximation (jlip . that uses 
the function and its first derivative to evaluate the fractional derivative, seems to be a good 
candidate in those cases. Suppose that we know the values of x(tj) onn + 1 distinct points in 
a given interval [a, b], i.e., for ti, i = 0, 1, . . . , n, with to = a and t n = b. According to formula 
([lip , the value of the fractional derivative of x(-) at each point ti is given approximately by 

N 

a D?x(ti) - A(a, N)(ti - a)- a x{U) + B(a, N)(t t - af^xiU) - £ C(p, a)(U - af^V^U) . 

p=2 

The values of x(ti), i = 0, 1, . . . , n, are given. A good approximation for x(ti) can be obtained 
using the forward, centered, or backward difference approximation of the first-order derivative 
[30] . For Vp{ti) one can either use the definition and compute the integral numerically, i.e., 
V p (ti) = r**(l — p)(t — a) p_2 x(r)dr, or it is possible to solve (fTUj) as an initial value problem. 
All required computations are straightforward and only need to be implemented with the 
desired accuracy. The only thing to take care is the way of choosing a good order, N, in 
the formula (jlip . Because no value of N, guaranteeing the error to be smaller than a certain 
preassigned number, is known a priori, we start with some prescribed value for iV and increase 
it step by step. In each step we compare, using an appropriate norm, the result with the one 
of previous step. For instance, one can use the Euclidean norm \\( a Df) new — { a Df) old \\2 
and terminate the procedure when it's value is smaller than a predefined e. For illustrative 
purposes, we compute the fractional derivatives of order a = 0.5 for tabular data extracted 
from x(t) = t and x(t) = e 2t . The results are given in Figure El 

4 Numerical solution to fractional differential equations 

The classical theory of ordinary differential equations is a well developed field with many tools 
available for numerical purposes. Using the approximations (|4|) and (jlip . one can transform 
a fractional ordinary differential equation into a classical ODE. 

We should mention here that, using @, derivatives of higher-order appear in the resulting 
ODE, while we only have a limited number of initial or boundary conditions available. In 
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Figure 6: Fractional derivatives of tabular data 



this case the value of N, the order of approximation, should be equal to the number of given 
conditions. If we choose a larger N, we will encounter lack of initial or boundary conditions. 
This problem is not present in the case in which we use the approximation (jllj) . because the 
initial values for the auxiliary variables V p , p = 2, 3, . . ., are known and we don't need any 
extra information. 

Consider, as an example, the following initial value problem: 



D°- 5 x{t) + x(t) = t 2 + j^ti, 
x(0) = 0. 



(20) 



We know that o-D?' 5 (i 2 ) = r ^ 5 ) ^ 2 • Therefore, the analytic solution for system ([20]) is x(t) = 
t 2 . Because only one initial condition is available, we can only expand the fractional derivative 
up to the first derivative in Q . One has 



1.5642 t-°- 5 x(t) + 0.5642 t°- 5 x(t) = t 2 + 1.5045 t 1 *, 
x(0) = 0. 



(21) 



This is a classical initial value problem and can be easily treated numerically. The solution 



is drawn in Figure 7(a) . As expected, the result is not satisfactory. Let us now use the 
approximation given by (|lll) . The system in (|20j) becomes 



A(N)t-°- 5 x(t) + B(N)t°- 5 x(t) - Zp =2 C(p)t°- 5 -PV p + x(t) =t 2 + r^fgyt 1 - 5 , 
V p (t) = (l-p)(t-a)P- 2 x(t), p = 2,3,...,N, 
x(0) = 0, 
[ V p (0) = 0, p = 2,3,...,N. 



(22) 



We solve this initial value problem for N = 7. The Matlab ode45 built-in function is used to 
integrate system (|22p . The solution is given in Figure [7(b)| and shows a better approximation 
when compared with (|21j) . 
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(a) Exact versus Approximation Q (b) Exact versus Approximation 



Figure 7: Two approximations applied to fractional differential equation (l20|) 



Remark 4.1. To show the difference caused by the appearance of the first derivative in 
formula (|lip . we solve the initial value problem (|20j) with B(a,N) = 0. Since the original 
fractional differential equation does not depend on integer order derivatives of function x(-), 
i.e., it has the form 

a D?x(t) + f(x,t) = 0, 

by (|14p the dependence to derivatives of x(-) vanishes. In this case one needs to apply the 
operator a Dl~ a to the above equation and obtain 

x(t) + a Dt~ a [f(x,t)]=0. 

Nevertheless, we can use (jlip directly without any trouble. Figure shows that at least for 
a moderate accurate method, like the Matlab routine ode45, taking B(a,N) ^ into account 
gives a better approximation. 



5 Application to the fractional calculus of variations 

The fractional calculus of variations consists in the study of dynamic optimization problems 
in which the objective functional or constraints depend on derivatives and/or integrals of 
non-integer order. This is a recent and promising research subject, under strong current 
research (see, e.g., [14 } 120 1 f2Tj and references therein). Here we show how to use expansions 
to transform a fractional problem into a classical one, where we can benefit from the vast 
number of techniques available in the field. Consider the following fractional variational 
problem including a left Riemann-Liouville fractional derivative, a Df , of order a G (0, 1): 

min J[x(-)] = J L(t,x(t),x(t), a D?x(t))dt 
x(a) = x a , x(b) = x b . 

One can deduce fractional necessary optimality equations to problem ([23]) of Euler-Lagrange 
type \22\ 123] : if x(-) is a solution to problem (]23p . then it satisfies the fractional Euler- 
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Figure 8: Comparison of our approach to that of [7J 

Lagrange equation 

dL „ dL d (dL\ . . 

There are a few attempts in the literature to present analytic solutions to fractional 
variational problems. Simple problems have been treated in [T]; some other examples are 
presented in [3] . In this work we use the two expansions discussed in Section [2] to reduce a 
fractional problem to a problem with derivatives of integer order. In order to illustrate the 
usefulness of our ideas in the area of the calculus of variations, we need to consider problems 
(123 p with a known exact solution. Examples 15.11 and 15.21 below are suitable for our purposes, 
since the analytic solutions can be easily obtained. Knowing the exact solutions, we compare 
the effectiveness of different approximation methods. 

Example 5.1. Let a £ (0, 1). Consider the following minimization problem: 



min J[x(-)}= [ [ Q D?x{t)-x 2 {t)\dt 
Jo 

s(0)=0, x(l) = 1. 
In this case the Euler-Lagrange equation (|24p gives 



(25) 



t Dfl + 2x(t) = 0, or x(t) = -— i -(1 - t)- a , 

IV (1 — a) 

which subject to the given boundary conditions has solution 

^ = -w^ 1 ~ t ^' + ( 1 --m^) t+ w^j- (26) 

The Lagrangian in Example 15.11 is linear with respect to the fractional derivative. This 
linearity makes the fractional Euler-Lagrange equation easy to solve. In a slightly different 
situation, e.g. Example 15.21 there is no well-known methods to solve the Euler-Lagrange 
equation (|2T 
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Example 5.2. Given a G (0, 1), consider now the functional 

J[x(-)} = f\ D?x(t) - Ifdt, (27) 
J o 

to be minimized subject to the boundary conditions x(0) = and x(l) = ^(a+i) ■ Since the 
integrand in (f27l) is non-negative, the functional attains its minimum when $Dfx{t) = 1, i.e., 

5.1 Numerical solutions to Example 15.11 

We use two different approaches. 

5.1.1 Expansion to integer orders 

Using approximation (jj]) for the fractional derivative in (|25p , we get the approximated problem 

N 



min J[a;(-)] = / 
Jo 



^C(n,a)t n - Q x (n) (t) -x 2 (t) 



_n=0 

x(0) = 0, x(l) = 1, 



dt 

(28) 



which is a classical higher-order problem of the calculus of variations that depends on deriva- 
tives up to order N. The corresponding necessary optimality condition is a well-known result. 

Theorem 5.3 (cf., e.g., [IB)). Suppose that x(-) 6 C 2N [a,b] minimizes 

rb 

L{t, x(t),x {1) (t),xM (t),..., (t))dt 
wi/i given boundary conditions 

x(a) = a , x(6) = 6 , 
x (1) (a) = ai, x (1) (6) = 6i, 

x^" 1 )(a) = a Ar _ 1 , x( JV - 1 )(6) = 6 JV _ 1 . 
T/ien x(-) satisfies the Euler-Lagrange equation 

dx dt \dx(v J dt 2 \dxW J { ' dt N \d x w ) 1 ' 

In general (j29[) is an ODE of order 2N, depending on the order N of the approximation 
we choose, and the method leaves 2N — 2 parameters unknown. In our example, however, 
the Lagrangian in (J28]) is linear with respect to all derivatives of order higher than two. The 
resulting Euler-Lagrange equation is the second order ODE 

N _ , 

^(-l)"C(n, «)^r(i n - a ) - I [-2i(t)] = 

n=0 
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Figure 9: Analytic vs. approximate solutions to Example 15.11 using approximation @. 



that has solution 



x(t) 



1 



2r(3 - a) 



N 



£(-l)T(n + l-a)C( 



n, a 



.n=0 



,2-a 



n, a 



n=0 



Figure [9] shows the analytic solution together with several approximations. It reveals that by 
increasing N, approximate solutions do not converge to the analytic one. The reason is the 
fact that the solution (|26p to Example 15.11 is not an analytic function. We conclude that Q 
may not be a good choice to approximate fractional variational problems. In contrast, as we 
shall see, the approximation (jlip introduced in this paper leads to good results. 



5.1.2 Expansion through the moments of a function 

If we use (jlip to approximate the optimization problem (|25p , we have 



J[x(-)] 



N 



A{a, N)t~ a x{t) + B(a, N^^x® - C{a,p)t x - p - a V p {t) - x 2 {t) 

p=2 

V p (t) = {l-p)€- 2 x{t), p = 2,3,...,N, 
V p (0) = 0, p = 2,3,...,iV, 
x(0) = 0, x(l) = l. 



dt, 



(30) 



Problem (|3U|) is constrained with a set of ordinary differential equations and is natural to look 
to it as an optimal control problem [25j. For that we introduce the control variable u(t) = x(t). 



17 




Figure 10: Analytic vs. approximate solutions to Example 15.11 using approximation (fTTj) . 



Then, using the Lagrange multipliers Ai, A2, . . . , Aj\r, and the Hamiltonian system, one can 
reduce (j30j) to the study of the two point boundary value problem 

±B(a, iV)* 1 -" - |Ai(t), 
(1 -p)tP- 2 x(t), p = 2,3,...,N, 



) 5 

2 < 



^(a,iV)t^-E; V =2 (l-p)^ 2 A P (i), 
-C(a,p)t( 1 -P- a ), p = 2,3, 



(31) 



AT, 



Ai(t) 
I Ap(t) 
with boundary conditions 
a?(0) = 0, 

^(0)=0, p = 2,3,...,JV, 

where cc(0) = and x(l) = 1 are given. We have V p (0) = 0, p = 2, 3, . . . , N, due to (fTU|) and 
A p (l) = 0, p = 2, 3, . . . , N, because V p is free at final time for p = 2, 3, . . . , N [25]. In general, 
the Hamiltonian system is a nonlinear, hard to solve, two point boundary value problem 
that needs special numerical methods. In this case, however, (]3ip is a non-coupled system of 
ordinary differential equations and is easily solved to give 



Ap(l) 



1. 



0, p = 2,3,... ,N, 



x(t) = M(a,N)t 



2-a 



N 

E 

p=2 



C(a,p) 



2p{2 — p — a 



■t p 



N 



C(a,p) 



p=2 



2p(2 - p 



a 



t, 



where 



M(a,N) 



1 



2(2 - a) 



B(a,N) 



A(a,N) 



N 



1 



O 



En 

p=2 



C(a,p)(l-p) 
(l-a)(2-p-a) 



Figure [10] shows the graph of x(-) for different values of N. 



5.2 Numerical solutions to Example 15.21 

Similarly to Example 15. 1[ it turns out that expansion 
while (fTT]) leads to good results. 



does not provide a good method 
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5.2.1 Expansion to integer orders 

Using dU) as an approximation for the fractional derivative in (|27|) gives 

2 



min J[x(-)] = (^2C{n,a)t n - a x (n) (t) - lj dt, 

s(0)=0, x(l)= 1 

T(a + 1) 

The Euler-Lagrange equation (I29p gives a 2N order ODE. For N > 2 this approach is 
inappropriate since the two given boundary conditions x(0) = and x(l) = r (a+i) are n °t 
enough to determine the 2N constants of integration. 

5.2.2 Expansion through the moments of a function 

Let us approximate Example 15.21 using ([lip . The resulting minimization problem has the 
following form: 



dt, 



2 

min J[x(-)] = / A(a, N)t- a x(t) + B(a, N^^x^) - V C(a,p)t l - p - a V p (t) - 1 
° I p =2 
V p (t) = (l-p)tP- 2 x(t), p = 2,3,...,N, 
V p (0) = 0, p = 2,3,...,N, 

x(0) = 0, x(l)= 1 

F[a + 1) 

(32) 

Following the classical optimal control approach of Pontryagin [25j as in Example 15.11 this 
time with 

N 

u(t) = A(a, N)t~ a x(t) + B(a, TV)* 1 ""^*) - ^ C(a,p)t 1_1, - a, ^,(t), 

p=2 



we conclude that the solution to ([32]) satisfies the system of differential equations 

' x(t) = —AB~H~ l x(t) + Zp=2 B^Cpt-PVpit) + iB-H^X^t) + B-H a -\ 
V p {t) =(l-p)tP- 2 x(t), p = 2,3,...,N, ' 

Ai(t) =AB-h- 1 \ 1 -Zp =2 (l-p)t p - 2 K(t), 
k \ p (t) =-B- 1 C(a,p)t-P\ 1 , p = 2,3,...,N, 

where A = A(a,N), B = B(a,N) and C p = C(a,p) are defined according to Section [21 
subject to the boundary conditions 

V p (0) = 0, p = 2,3,...,N, \A P (1) = 0, p = 2,3,...,N. { ' 

The solution to system (|33p - ([34p . with ./V = 2, is shown in Figure PTT1 
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Figure 11: Analytic versus approximate solution to Example 15.21 using approximation (fTT|) . 



6 Conclusion 

During the last three decades, several numerical methods have been developed in the field of 
fractional calculus. Some of their advantages, disadvantages, and improvements, are given in 
[3]. Based on two continuous expansion formulas (|2j) and (J9j) for the left Riemann-Liouville 
fractional derivative, we studied two approximations (jU) and (jlip and their applications in the 
computation of fractional derivatives. Despite the fact that the approximation @ encoun- 
ters some difficulties from the presence of higher-order derivatives, it exhibits better results 
regarding convergence. Approximation (jlip can also be generalized to include higher-order 
derivatives in the form of (I15p . The possibility of using () 1 If) to compute fractional derivatives 
for a set of tabular data was discussed. Fractional differential equations are also treated suc- 
cessfully. In this case the lack of initial conditions makes (j4| less useful. In contrast, one can 
freely increase N, the order of approximation (jlip . and find better approximations. Compar- 
ing with (PT4|) . our modification provides better results. We finished by discussing the solution 
of fractional variational problems using the introduced approximations. Similar methods can 
also be applied to fractional optimal control problems. 

For fractional variational problems, the proposed expansions may be used at two different 
stages during the solution procedure. The first approach, the one considered in Sections 15.11 
and 15. 2\ consists in a direct approximation of the problem, and then treating it as a classical 
problem, using standard methods to solve it. The second approach would be to apply the 
fractional Euler-Lagrange equation and then to use the approximations in order to obtain 
a classical differential equation. However, the Euler-Lagrange equations (|24p involve right 
Riemann-Liouville derivatives, which introduces undesirable difficulties. 
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